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Abstract. Results from helically symmetric scalar field models and first results 
from a convergent helically symmetric binary neutron star code are reported here; 
these are models stationary in the rotating frame of a source with constant angular 
velocity Q. In the scalar field models and the neutron star code, helical symmetry 
leads to a system of mixed elliptic-hyperbolic character. The scalar field models 
involve nonlinear terms of the form ip 3 , (Vi/>) 2 , and if>Oip that mimic nonlinear 
terms of the Einstein equation. Convergence is strikingly different for different 
signs of each nonlinear term; it is typically insensitive to the iterative method 
used; and it improves with an outer boundary in the near zone. In the neutron 
star code, one has no control on the sign of the source, and convergence has 
been achieved only for an outer boundary less than ~ 1 wavelength from the 
source or for a code that imposes helical symmetry only inside a near zone of that 
size. The inaccuracy of helically symmetric solutions with appropriate boundary 
conditions should be comparable to the inaccuracy of a waveless formalism that 
neglects gravitational waves; and the (near zone) solutions we obtain for waveless 
and helically symmetric BNS codes with the same boundary conditions nearly 
coincide. 
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1. Introduction 

Initial data for the inspiral of binary neutron star (BNS) systems and corresponding 
quasiequilibrium BNS models have been based either on the initial value equations 
alone or on the IWM (Isenberg- Wilson-Mathews) spatially conformally flat ansatz. In 
each case one solves a truncated version of the Einstein equation for a metric having 
fewer than the six independent potentials that remain after a choice of gauge. The 
error of the approximation limits the accuracy of the first orbits of simulations in two 
ways: The initial data has spurious radiation; and, more importantly, the balance 
between gravitational attraction and orbital acceleration is not enforced, leading to 
orbits that are not exactly circular. 

One way to go beyond spatial conformal flatness is to construct an analog in 
the full theory of the Newtonian binaries that are stationary in a rotating frame. In 
general relativity, these models are helically symmetric spacetimes 0GJi with equal 
amounts of ingoing and outgoing radiation. Binary black holes of this kind were first 
discussed by Blackburn and Detweiler j^j, and models involving nonlinear scalar wave 
equations have been studied by a group of researchers organized by Price (henceforth 
the Consortium) H El EIEI El El Because the power radiated by a helically symmetric 
binary is constant in time, the spacetime cannot be asymptotically flat. At distances 
large compared to l/f2, however, the spacetime of a binary system with a helical Killing 
vector approximates asymptotic flatness - until, beyond about 10 M for neutron-star 
models of mass M, the enclosed energy in gravitational waves is comparable to the 
mass of the binary. The approximation is similar to, and possibly more accurate than, 
a 3rd post-Newtonian approximation in which the 2 1/2 post-Newtonian radiation is 
omitted. 

In this paper, we report the construction of a convergent neutron-star code in 
which the full Einstein equation, together with the equation of hydrostatic equilibrium, 
is solved numerically under the assumption of helical symmetry. Convergence relies 
on a boundary that is not much farther than a wavelength from the system, and 
we present results from a number of related nonlinear scalar field models in which 
convergence requires either a small coefficient of the nonlinear term or a boundary 
close to the source. The results of these toy models are surprising in two ways. First, 
convergence of the scalar-field models is most strongly affected by the sign of the source 
term, with one choice of sign yielding a convergent solution for remarkably large values 
of the nonlinear terms we examined. Second, convergence does not depend strongly 
on the iterative method used to solve the equation, on whether one uses, for example, 
a Newton-Raphson iteration or an iteration based on a Green function that inverts 
only a convenient part of the second-order nonlinear operator. 

The plan of the paper is as follows. Sect. 2 introduces the set of toy scalar field 
models with nonlinear terms chosen to mimic the nonlinear terms in the dynamical 
part of the Einstein equation. We describe several iteration methods, one closely 
related to that of the Uryu code, the others to codes developed by Andrade et. al. 
in \J\ and further by Bromley, Owen, and Price HJ. Sect. compares the accuracy 
of codes based on the various iteration methods. Sect. 0] reports the main results on 
convergence of the scalar-field models. Finally, Sect. El describes Uryu's neutron-star 
code, presents its first results, and compares the solution it yields to that of a closely 
related code based on a waveless formalism ^OJ ■ 
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2. Toy Problems 

2.1. Helically symmetric, nonlinear scalar wave equations 

We consider a scalar field tp on Minkowski space, satisfying a wave equation with a 
source s that mimics two objects in circular orbit and with a nonlinear term Af[ip], 
whose strength is adjusted by a coefficient A: 

Dip - AA/ty] = s. (1) 

We use three different nonlinear terms: JV[ip] — V' 3 , A/"[i/'] = |Vi/'| 2 with V the spatial 
gradient, and Af[ip] = ip^^i chosen to represent the types of nonlinear terms that 
appear in dynamical components of the field equations. 

The source s is a sum of two 3-dimensional gaussian distributions, 




centered about points ±R, R(£) = a [cos(f2t)x + sin(£!t)y], each a distance a from the 
origin and each having spread a 2 and total charge q. The source s is stationary in a 
frame moving with angular velocity f2; that is, it is Lie-derived by the helical Killing 
vector 

k a = t a + n<j) a (3) 

of Minkowski space, where t a and <p a (equivalently dt and d^) are generators of time- 
translations and of rotations in the plane of the binary source. 

One can regard the scalar-field models as toy models of neutron stars of mass M, 
if the charge q of each gaussian source is identified with 4irM. In gravitational units 
(G=c=l), all quantities of a binary star system can be specified in terms of M. In 
the models presented below, we set q = 1, a = 1, a = 0.5, and fl = 0.3, corresponding 
to a binary system of mass M, stellar separation 2a = 8ttM, stellar radius a — 2ttM, 
and velocity v — afl = 0.3. 

A helically symmetric solution, like a genuinely stationary solution, is given 
by its value on a spacelike slice and the field equation can be written in a 
form that involves only spatial derivatives. That is, using the symmetry relation 
£fc"0 = {dt + fld^'ip — to replace time derivatives by <f> derivatives, one can rewrite 
Eq. at t = in the spatial form, 

(v 2 - n 2 dl)t> - W[*] = s, (4) 

where ^(r,6»,0) = ip(t = 0, r, (9, </>), S{r,6,4>) = s{t = O,r,6»,0). Then ij){t,r,6,4>) = 
#(r,0,0-Qi). 

The operator C := V 2 — ^ 2 9^ has a mixed character, elliptic inside the light 

cylinder Q^J x 2 + y 2 = 1, hyperbolic outside. The difficulties in finding numerical 
solutions stem from this behavior. In finding an iterative solution, one inverts the 
operator C, but C is not negative, and it lacks the contraction property that underlies 
the convergence of iterative schemes used to invert nonlinear elliptic equations and to 
prove existence of exact solutions. 
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2.2. Numerical methods 
KEH method 

In the KEH method one splits off the linear, flat-space operator C and inverts it 
by a choice £ _1 of Green function. The iterative solution, beginning with "Jo = C _1 S, 
is then given by 

*n+l = *0 + AT-W[* B ] (5) 

Although C is not elliptic, the operator associated with each spherical harmonic 

is the elliptic Helmholtz operator: {r)Y lm \ = [V 2 + m 2 Cl 2 ][^i m (r)Yi m ]. The 

d 2 2 d 1(1 + 1) 2n2 , . . , 

corresponding radial operator — r- H - \-m 11 nas as its eigentunctions 

dr 1 r dr r 1 

the spherical Bessel functions, from which a Green function is constructed. For 
defmiteness, we choose as £ _1 the form that, for a bounded source, yields the half- 
advanced+half-retarded solution, namely 

C^S :=5Z/ dr'S lm {r')g lm (r,r')Y lm (0^), (6) 



with 



1 rK 

910 = o; ; : -7XD gi m = mQ,ji(mVLr < )ni(mQ,r > ), m/0. (7) 
2t + 1 r^J 

Here r< = min(r, r'), r> = max(r,r'), and jz(x) and n;(x) are the spherical Bessel 
functions of the first and second kinds. 

At each iteration, the nonlinear term A/"[vI/ n ] serves as an effective source. The 
polynomial nonlinear function \1/ 3 is most easily computed by shifting from ^i m (r) 
back into \&(r, 9, <fi) and cubing at each point, while the derivative-based nonlinear 
terms, IV^I 2 and are calculated using the properties of the spherical harmonics. 

Finally, as is usual in codes to solve nonlinear elliptic equations, we use softening 
and continuation to extend the range of convergence to larger values of A. That is, 
instead of using ^n+i as defined in Eq. 10), we can use a softened ^^+1 defined by 

= w^ n+x + (1 - w)tf n . (8) 

Given a converged field solution for some nonlinearity with small A, it is sometimes 
possible to use continuation to obtain a solution for larger A: The converged solution 
to Eq. with weak nonlinearity is used as the initial field \&o for the iteration of 
Eq. Q for strong nonlinearity. In this way one moves along a sequence of solutions 
with increasing values of A. The effectiveness of softening and convergence is explored 
in Section 01 

Finite difference and eigenspectral methods 

The Finite Difference code uses an iteration based on the Newton-Raphson 
method, with numerical approximations that reduce it in part to a secant method. 
Write the equation to be solved as 

F = Cf- -5 = 0. (9) 

Numerically, \I/ is given by a set of values on the three-dimensional spatial slice. 
Given an initial field 'J,;, each iterative step generates a modification S^i by inverting 

Jii&Vj = -Ft (10) 

where is the Jacobian 
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The Helmholtz operator has the form (£*f?)i = iy^j where is constructed 
from finite difference operations and incorporates boundary conditions. The 
corresponding part of Jy is simple. 



Jij = [£«*k - A(JV[*])< - S«] (12) 

A ~9*7" (13) 

The nonlinear piece of the Jacobian is evaluated numerically by varying local field 
values. 

The Eigenspectral code J0J uses the same iterative scheme as the Finite Difference 
method, but it employs adapted coordinates and a discretized spectral decomposition. 
To specify the adapted coordinate system, we begin with comoving Cartesian 
coordinates (x,y,z) where the z-axis is the axis of rotation. The axes are rotated 
to a set 

X = y,Y = z,Z = x (14) 

for which the Z-axis is a line through the center of each source. The adapted 
coordinates are chosen to approach spherical polar coordinates, with measured 
from the Z axis, far from the sources. For each point (X, Y, Z) , let r% and r-i be the 
distances from the source centers, 0\ and 6*2 corresponding angles from the Z axis. 
Then, 



n = J{Z -a) 2 + X 2 + Y 2 , (15) 



r 2 = \J{Z + a) 2 + X 2 + Y 2 . (16) 
The adapted coordinates are 

X = V r i r 2 (17) 

$ = tan-^X/Y) (19) 
The spectral decomposition involves the angular Laplacian of and $, 



L> 2 $ := 



_L d_ 

sin 90 



• n 9 



1 d 2 

(sin0) 2 9$2' 



(20) 



note that D 2 is not the angular part of V 2 in adapted coordinates, but it 
agrees asymptotically with the usual angular Laplacian far from the source. 
Instead of the spherical harmonics of the continuum Laplacian D 2 , the 
eigenspectral code uses the exact eigenvectors of the matrix L obtained by 
discretizing D 2 on the adapted coordinate grid {0j,$j}. Angular derivatives 
are represented in L by second order finite differencing. That is, with 
[sin0£> 2 $] (e o ,$ 6 ) approximated by £V ■ L a b,ijV(®i, 

(k) 

the normalized eigenvectors of L satisfy 

E^M^f -AWataerf, (21) 
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(22) 



With the field expanded in terms of the eigenvectors Y^- 



*(x,e <) * i ) = X;a ( * ) (x)^ 



(23) 



k 



and the operator £ written in terms of adapted coordinates (Eq. 8-17 of 0), >C\& is 
expressed in terms of the Yy (Eq. 31 of |5]), 



where the ay Pk' ,k, and Jk',k involve angular derivatives of the Y^- 1 computed by 
finite differencing (Eqs. 41-43 of 9 ) . This equation for is iterated in the same 
fashion as in the Finite Difference method to find the solution with nonlinear terms. 

If all harmonics were retained, the Eigenspectral method would be the equivalent 
of the Finite Difference method in adapted coordinates, although in a different 
basis. The advantage of the adapted coordinates is that the distribution of points 
encodes most of the physically relevant information in the low-order harmonics. Its 
disadvantage is that higher order harmonics require an increasingly cumbersome 
formalism. The code is consequently limited to harmonics I < 2, saving enough 
memory to allow high resolution in the radial coordinate. 

2.3. Boundary conditions 

The KEH method at each iteration finds a solution that has standing wave behaviour 
for the flat space piece, imposing boundary conditions by the choice of Green function. 
There remains freedom to add a homogeneous solution at each iteration, and this has 
been used to study the sensitivity of convergence to the choice of boundary condition. 

The Eigenspectral and Finite Difference methods include boundary conditions in 
the finite difference matrix for the linear □ operator. Outgoing conditions are imposed 
on the edges of the grid by enforcing 



A solution for ingoing radiation can then be generated by a spatial inversion across 
the plane through the sources and perpendicular to their rotation. At each step, a 
periodic solution is constructed by superposition of the ingoing and outgoing solutions. 

3. Estimating numerical accuracy 

3.1. KEH code 

Several types of numerical approximation in the KEH code produce effects that can be 
estimated by convergence testing. Most obvious is the choice of spatial grid in (r, 6, <p) 
on which numerical integration is performed. Very high resolution in 9 and <f> is easily 
obtained. The number of radial points is more problematic with our simple equispaced 
grids. A Richardson extrapolation error estimate from varying radial grid spacing 
shows that precision of 10 -5 is estimated for runs comparing code results. Estimation 
of the range of A giving convergent solutions is done at lower radial precision. 




(24) 



(d r V - 9<^)r=rw = 



(25) 
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Figure 1. The rms change in the converged field as the number of harmonics 
considered increases. 



As the KEH method rests on the decomposition of the field into spherical 
harmonics, accuracy will depend on the number of harmonics retained in the numerical 
calculation. Examining the difference between results with an increasing number of 
harmonics retained shows that a good approximation is to use up to I = 12 in the 
code. The fractional difference between the field calculated with ^ max = 12 and the 
field including higher harmonics is less than 10~ 6 at each point. Figure ^ displays 
the rms difference in solutions as the number of harmonics retained by the code is 
increased. 

In the toy models presented below, a solution is computed at each iteration using 
the half-advanced+half-retarded Green function 10. Because the linear field $ of a 
perpetually radiating source falls off like r -1 , the nonlinear terms TV = (V^) 2 and 
J\f = that serve as effective sources for each iteration do not fall off fast enough 

for the integrals to converge, if the outer boundary extends to infinity. One can, 
however, pick out a solution that is independent of the outer boundary by fixing the 
value of at a finite radius R. That is, one can, at each iteration, add the homogeneous 
solution that maintains the specified value of \& at R. This has remarkably little 
effect on the field in the region close to the sources, with less that a 1% change 
in field strength for points with r < 6. This insensitivity of the near-zone field to 
the amplitude of the waves, when the source dominates the solution, is the reason a 
helically symmetric solution makes sense as an approximation to an outgoing solution. 
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3.2. Comparing codes 

With different 3D grid patterns for the codes, it is most straightforward to compare 
results on rays through the volume of interest. We compare the results extrapolated to 
three orthogonal axes: taking the x-axis through the centers of the two sources and the 
z-axis perpendicular to the plane of rotation. Both x and y axes show wave behaviour 
away from the sources; along the z axis, tp shows only Coulomb-type behavior, because 
Yim vanishes on the axis when m ^ 0, and for m = 0, Cip = V 2, 0- 

Sample comparison plots are shown in Figure [3 Further comparisons are shown 
in the longer version of this paper on the gr-qc archive. 

The field values on the axes are interpolated for comparison. Differences between 
fields are divided by the average field value of the three codes to find relative error, 
as plotted in the insets of Figure [21 We compute the rms of this relative error for the 
grid points along each axis. These rms values are expressed as percentages in Table 
^ In the worst case, the rms difference is 3% between codes. 

A smaller rms difference on the rotation axis than on the source and perpendicular 
axes indicates that discrepancies in the wave region dominate, as in the FD-KEH 
comparison with Af[^f] = , 5 3 , A = 100. In other cases, the error on all three axes is 
comparable, and there is some shift between codes seen even on the rotation axis, as 
in the FD-KEH comparison with Af[$] = |V*| 2 , A = 3. 

Table 1. A comparison of code output for scalar models. The values in the table 
give the percent rms difference between of if-values from three numerical methods 
as in Figure 1. The rms differences were computed from the points along each 
principal axis. 





Source axis 


Perp. axis 


Rotation axis 


Linear (jV[tf] = 0) 








FD to KEH 


0.25 


0.24 


0.10 


FD to ES 


1.02 


1.09 


1.08 


ES to KEH 


1.07 


1.15 


1.11 


= * 3 , A = 100 








FD to KEH 


3.56 


3.44 


0.21 


FD to ES 


1.22 


1.67 


1.14 


ES to KEH 


3.28 


2.86 


1.19 


Affif] = f 3 , A = -2.2 








FD to KEH 


0.35 


0.36 


0.21 


FD to ES 


1.96 


2.02 


2.08 


ES to KEH 


1.77 


1.99 


1.92 


M\9] = |W| 2 , A = 3 








FD to KEH 


1.73 


1.47 


1.33 


FD to ES 


1.39 


1.80 


1.44 


ES to KEH 


2.69 


2.76 


2.77 


7V[*] = |W| 2 , A = -100 








ES to KEH 


1.96 


2.20 


1.65 


JV[*] = A = -0.7 








FD to KEH 


0.53 


0.48 


0.20 


FD to ES 


1.06 


1.13 


1.13 


ES to KEH 


0.99 


1.23 


1.00 
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Figure 2. The scalar field <t as a function of distance r from the origin in units 
of the orbital radii along the source axis. Each panel corresponds to a model 
specified by the form of the nonlinear term AT[*i?] and the value of the parameter 
A, as written above the panels. In all cases, the angular frequency of rotation is 
f! = 0.3 and the source strengths are unity. The plots show results from the KEH 
method (solid red curves), the Eigenspectral method (dashed green) and the Finite 
Difference code (dotted blue). The insets give more detailed comparisons between 
these results; the curves show the difference between pairs of solutions relative to 
the average field. The Eigenspectral method gives a stronger periodic modulation 
relative to the other methods, due to its restriction to low-order harmonics. 
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Table 2. The range of nonlinear amplitude A in scalar models for which codes 
converged. The first column gives the type of nonlincarity, while the second 
column indicates the numerical method, as discussed in the text. The third and 
fourth columns, showing A m i n and A max , give the range of A for which convergence 
was achieved. Where an inequality is given, no limiting value was found. 





Code 


•\nin 


^max 


$3 


KEH 


-2.3 


n 


$3 


KEH, softened/continued 


-2.3 


>10000 


$3 


FD 


-2.5 


>1000 


$3 


ES 


-2.4 


>1000 


|W| 2 


KEH 


-30 


7.7 


|V*| 2 


KEH, softened/continued 


<-1000 


7.7 


w 2 


FD 


-2.0 


5.6 


w 2 


ES 


-360 


7.3 


w 2 


ES, continued 


<-1000 


7.3 




KEH 


-0.96 


1.3 




KEH, softened/continued 


-0.96 


10 




FD 


-1.8 


>1000 




ES 


-1.7 


>1000 



4. Measuring range of convergence 

For each code, a series of runs varying the nonlinear amplitude A determines the 
range of A for which the code gives a convergent solution. The source distribution and 
boundary location are fixed. 

Softening and continuation are used to extend the maximum absolute value of 
A that allows convergence. For each nonlinear term, strikingly different behaviour is 
found for opposite signs of A. That is, denoting by A max and A m i n the largest positive 
value and the smallest negative value of A for which a code converges, we find for each 
nonlinear term that A max and |A m i n | differ by at least a factor of 100. In each case, 
the sign of A that yields greatest convergence is opposite to the sign of the source 
term, where the source is large. The 'favorable' sign damps the effect of the source 
distribution on the field, while the 'unfavorable' sign gives an amplified effective source. 
For 7V[*] = * 3 , |V*| 2 ,*D*, the favorable signs of lambda are +,—,+, respectively. 

Softening and continuation improved the range of convergence for the favorable 
sign, but had no effect on the limit of the unfavorably signed A. Figure [3] shows 
the effect of softening, parametrized by lo, and of using continuation on the limiting 
favorable-sign values of A. 

Convergence is attained, when softening and continuation are used, for all 
attempted favorable sign values for A, except when using the KEH method with the 
^CM nonlinear term or when using the Finite Difference method with the | X7\I> | 2 
nonlinear term. 

The range of convergence results are given in Table |21 If softening and 
continuation maintain convergence to the largest tested value of A, the uncertainty 
in the true maximum value or minimum value of A is indicated by > or <. 
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(b) 



rcontinuation 
continuation 
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no continuation 
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Figure 3. The effect of a softening parameter ui and continuation on range of A. 
The plots indicate the limiting value of A for convergence of the KEH code. Each 
panel corresponds to a different nonlinear model, with A/'^] set to Nl" 3 (panel a), 
IV^I 2 (panel b), and ^CM/ (panel c). In each case, the limiting A value of the 
opposite sign was unaffected by softening or coninuation. 



Boundary and convergence 

As described in Sect. below, Uryu has obtained a convergent helically symmetric 
BNS code within a region that extends about one wavelength from the center of the 
source, by using a waveless formulation outside that radius. Led by this result, we 
explored the effects of the outer boundary placement on the range of converging A for 
the fully helical code. 

While bringing the boundary in has little effect on the limits of the favorably 
signed A, it does allow a significantly larger magnitude of the unfavorable A, as shown 
in Fig. H 



5. Hybrid helical/ waveless equations 

One expects a helically symmetric solution "J to accurately approximate the outgoing 
solution only in the near zone, where the Coulomb part of the field is dominant. 
That is, the two solutions agree when one can ignore the wavelike character of ^. 
In that situation, one should be able to approximate the outgoing solution with 
comparable accuracy by using a waveless approximation in which only the Coulomb 
part of the field is retained. We can obtain an approximation of this kind by replacing 
the D'Alembertian by the Laplacian - by discarding second time derivatives. In 
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Figure 4. The effect of the outer boundary R on convergence in the cubic 
nonlinear model (A/1 1 !/] = \E" 3 ). The limiting value of A for which the KEH method 
converges as a function of outer boundary location, R, in units of orbital radii. 



the last part l|5.1J) of this Section, we compare results of this version of a waveless 
approximation to a solution that has exact helical symmetry. We impose equivalent 
boundary conditions on the solutions by imposing helical symmetry only in the near 
zone, in effect matching to an exterior waveless solution. As we describe Sect. in 
Uryu's neutron-star code, a helically symmetric near-zone is similarly matched on an 
initial hypersurface to an exterior waveless solution to the Einstcin-Eulcr equations. 

We find as follows a helically symmetric near-zone VP matched to an exterior 
waveless solution. As before, after imposing helical symmetry, we have a scalar 
equation for the toy scalar model in the form 



where the effective source is a sum of the inhomogcncous source term and a nonlinear 
source. Instead of inverting the operator on the left hand side in the equation with 
the wavelike Green's function Cr 1 of Eq. JBJ, we split the operator, writing 



We regard the extra term on the right hand side as a part of the effective source, and 
invert the Laplacian operator to solve iteratively for ^ . 

We can now truncate the extra source term £l 2 d 2, $> (as in Uryu et al.^HI) a t 
some radius -Rtrunc written in units of the wavelength n/fl of the quadrupole wave. 
The resulting solution is helically symmetric inside the truncation radius and waveless 
outside. 

If this truncation radius is too large, the iteration fails to converge. In fact, this 
iterative method, in which we invert the Laplacian instead of the full second-order 
operator, fails to converge even for nonlinear terms that yield convergent solutions 
when we impose helical symmetry everywhere and use the iterative methods of the 
previous section. For a given small truncation radius, the critical value of A depends 
on the outer boundary radius i? ut, as in the fully helical result of Figure 0] 

In Fig. |5j we plot the critical value |A m i n | of A as a function of -Rtrunc, for the 
nonlinear term AvE' 3 with the unfavorable sign choice for A. As in the models discussed 
above, the inhomogeneous source term is the double Gaussian, with Q = 0.3. In 
the upper panel, we plot results of iteration without continuation. Different curves 
correspond to different value of outer boundary radius, i? ou t- As is clear from the 
figure, the value A cr = A m i n is nearly independent of -Rtrunc up to maximum value of 



(V 2 -n 2 dl)v = S + AA/M 



(26) 



V 2 * = S + AAf[¥] + n 2 d 2 ^. 



(27) 
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Figure 5. Critical parameter A m j n | is plotted as a function of iJtrunc- We plot 
curves for different computational outer boundaries -Rout- Note that iitrunc is 
normalized by 7r/Q. The range of convergence is given for the iterative method 
described in this section with continuation used to reach a limiting value of |A m j n |. 
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Figure 6. Left panel: Plot of hybrid hyb) and waveless (^wl) solutions 
along the source axis normalized by tt/Q, the wavelength of the I = m = 2 mode. 
Right panel: The absolute error between hybrid and waveless solutions along the 
same axis. 



-Rtrunc- Beyond that value, the iteration based on inverting the Laplacian does not 
converge for any finite value of A. This transition from convergence to non-convergence 
is abrupt, and we see a clear limit, -Rtrunc ~ 0.47, for the truncation radius. We also 
see that using continuation we reached an absolute limit of i?trunc ~ 0.47 for any value 

Of -Rout ■ 

We note that the truncation radius of the extra source (that is, Vl 2 d^S>) cannot 
be larger than a half of wave length of quadrupole component of wave. This is roughly 
consistent with the results of Uryu's code described below. 

5.1. Comparing hybrid and waveless results 

In Fig. I5.1l we plot the difference between a hybrid solution and a solution that is fully 
waveless. The profile is along the source (along the x-axis). The resulting similarity of 
the solutions confirms our expectation that waveless and helically symmetric solutions 
should agree in the near zone. The level of disagreement, between 0.05 and 2%, for 
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r less than one wavelength, appears primarily to result from the fact that the two 
solutions do not have identical boundary values at the boundary -Rtrunc 01 the region 
where one solution is helically symmetric and the other is waveless. That is, by 
truncating the d\ term in the effective source at r = i?trunc, we match the helically 
symmetric solution to an exterior waveless solution, but to a slightly different waveless 
solution than that obtained by using the waveless equation everywhere. 



6. Helically symmetric binary neutron star code 

In this section, we report the construction of a first helically symmetric code for 
binary neutron stars. The code modifies a recently developed numerical method for 
computing models of compact binary stars in a waveless approximation |1(J| to produce 
solutions to the full Einstein-Euler system (the Einstein-perfect fluid equations)with 
exact helical symmetry. Waveless and helically symmetric solutions are expected to be 
accurate in the near zone, where the gravitational wave amplitude is small compared 
to the Coulomb contribution to each metric potential; and in its present form, the 
helically symmetric code converges only in the near zone. Convergence of the code is 
achieved by using the waveless formulation outside one wavelength. In effect, the outer 
waveless solution imposes boundary conditions on the helically symmetric solution at 
the edge of the near zone. 



6.1. Formulation 

Our formulation is applied to Einstein's equation written in 3 + 1 form, with spacetimc 
M = K x S. Let n a be the future-pointing unit vector normal to the slices {t} x S 
of this foliation. The generator of the time coordinate t a and the helical symmetry 
vector k a — t a + Vt<t> a are related to n a and the shift (3 a in the usual way, 

t a = an" + P a and k a = an a + w Q , (28) 

where uj a is a rotating shift vector u a = (3 a + ri(j) a . 

The projection 4-tensor 7 Q ^ orthogonal to n a , defined by 

7q/3 = 9af3 + n a np, (29) 

is associated with the 3- metrics j a b(t) on the spatial slices St. (Any spatial tensor on 
S t can be identified with a spacetime tensor orthogonal to n a on all its indices). The 
metric then has the form, 

ds 2 = -a 2 dt 2 + jij (dx l + I3 l dt)(dx 3 + f3 j dt). (30) 

The extrinsic curvature of each slice E is defined by 

K a b = —-£ n 7ab, (31) 

where £ n operating on spatial tensors such as "f a b has the meaning 

£«7afc := —dtjab - -£/37afc, (32) 
a a 

with dt^ab the pullback of £t7a/3 to S. We introduce a conformal decomposition of 
the spatial metric, j a b = ip^Jab, and a condition 7 = /, where / is the determinant 
of the flat metric f a b and 7 the determinant of the conformal metric "/ a b- Further h a b 
and h ab are introduced by 

lab - fab + Kb , r b = f ab + h ab . (33) 
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Helical symmetry, £ k g a p = 0, implies for the 3-metric and extrinsic curvature, 
£ klab = 0, £ k K ab = 0. (34) 
Using the relation (|28|l between n a and k a , we have 

£n7afc = £^7ab> aim £nK a b = £u]K ab . (35) 

a a 
Projections of the Einstein equation along the normal n a are the Hamiltonian 
and momentum constraints, 

(G Q/3 - SnT a p)n a nP = 0, (36) 

(G al3 - 87^)7^ = 0. (37) 
while the spatial projection has trace and tracefree part 

(G aP - 8^T tt/3 ) 7 Q/3 = 0, (38) 

(G Q/3 - 8ttT^) ( 7a a 7b P - \ lab l a ^j = 0. (39) 

The above set of equations are solved for the conformal factor ip, the shift f3 a , 
the lapse a and the deviation of the conformal metric from the flat metric h ab , 
respectively, imposing coordinate conditions, the maximal slicing condition, K = 0, 

o 

and the generalized Dirac condition D b "f — 0, as in the waveless formulation [TT]. In 
the following section, we concentrate on the treatment of the spatially tracefree part 
of Einstein equation l|39l) . solved for h ab , in which the second time derivatives of h ab 
change the equation from elliptic (in the waveless formulation) to the mixed form that 
characterizes helically symmetric wave equations. A concrete treatment of the other 
components of the field equations that yield elliptic equations for ip, fl a , and a, as 
well as equations of matter source is given in Ref. 

The waveless code differs from the helically symmetric code by the requirement 
that jab have vanishing derivative along t a instead of vanishing derivative along 
k a . The extrinsic curvature and matter variables are still required to be helically 
symmetric. These requirements result in elliptic equations for the field variables, 
including the non-conformal part of the 3-metric h ab = j ab — f ab . 

In the helically symmetric code, with dt^ab nonzero, the elliptic equations are 
replaced by Helmholtz equations for h ab , whose source is almost identical to that of 
the waveless formulation. We implement a KEH solver for the helical formulation and 
investigate convergence of the iteration for a compact binary star source. 

6. 2. Helmholtz equation for h ab 
Eq. has the form 

(G a/3 - ^T a p) - = Sab - \labl Cd £cd = 0, (40) 

where 

E ab := -£ n K ab + 3 R ab + KK ab - 2K ac K b c - -D a D b a - &nS ab , (41) 

a 

with 3 R ab the Ricci tensor on S associated with 7 a b, and S ab the projection of the 
energy stress tensor, S ab := T af37a a ^ . 
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By isolating the terms, Oh ab := (— <9| + D c D c )h ab , that occur in £ ab in Eq. IgTJ, 
one can rewrite Eq. i|4U|) in the form 

□/lab = S ab , (42) 

o o 

where □ is the flat d'Alembertian operator and D a :— f D b . Then, as in the scalar 
models, helical symmetry of the conformal metric, £klab — £>kh a b — 0, results in the 
operator 

Dh ab = -0 2 t h ab + D c D c h ab = (A - n 2 £l)h ab , (43) 

o o o 

where the flat Laplacian is defined by A = D C D C . 

o 99 

Even when one uses the Cartesian components of h ab , however, A — f2 £^ does 
not coincide with the Hclmholtz operator, because d^h^ is a Cartesian component 

o o 

of <f> ' Dft, a fc, not of Lf/jhab- To isolate the Hclmholtz operator C = A — Q 9? = 
A — f2-(</> ■ D) , we write := <fr • D and find 

£^7a6 = <^7 afc + (<9 7 ac + £ 07ac ) i3 fc c + (<9 07cfc + £^ cb ) D a <j> c , (44) 

o o 

where a relation <jf D C (D a (jr) = is used. Moving all terms in Eq. I|44[l except d 2; y ab 
from the LHS to the RHS of Eq. |4"2")l. we obtain the Hclmholtz form 

Ch ab — S a b '■= 2 ^£ab ~ 7pabl Cd £cd 

- ^ ab D e h cd D e h cd + ^i a b^ 2 d^h cd d^h cd , (45) 
where the barred expression £ ab is defined by 

£ ab := iC + 3 Rt b - ^D a D b a ~ 2^ 4 A ac A b c - 8nS ab 

i A 



^ 9 

2a 



2 



(dflac + Efjac) D b (j> C + {d^ cb + E^Jch) D a (j) C 



if; 4 v 4 ~ ■0 s 

+ £/3£/37ah H A ab £ w ln — . (46) 

2ar a a 

o 

In this expression for £ ab , the coordinate conditions K = and D b ^f = 0, 
mentioned above, are imposed; and the tracefree part A ao of the extrinsic curvature, 
A ab := K ab — ^jahK, is introduced in the rescaled form A ab := ip~ 4 A ab . Terms R^ b 

and 3 Rf h arise from the conformal decomposition of the Ricci tensor |llj . 

The source (|46|l can be written concisely without separating the second derivative 
term &ih ab explicitly as above. Applying the helical symmetry condition (|35[) to 
Eq. 1411) and subtracting d^hab term from both side of Eq. Q40[). the source term of 
the exactly helical equation Eq. I|45|l can be rewritten 

£ ab -.= -£U^A ab )-l-n 2 dlh ab , 

a 2 v 

+ R^t + 3 Rt b - \p a D b a - 2i, A A ac A b c ~ ^S ab (47) 
which is equivalent to the above source term l|46[) . 
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6. 3. A different formulation for the use of elliptic solver, and a comparison with the 
waveless approximation 

Instead of isolating the Hclmholtz operator on the LHS for the case of helical 
symmetry, one can formally isolate an elliptic operator and solve the equation 
iteratively using an elliptic solver as in the waveless approximation. With this grouping 
of terms, the helical symmetry condition Ij35(l . applied to Eq. (|4 If) - leaves the term 
£u)K a b as part of the effective source, and the Laplacian of h a b is separated out from 

o 

3 R a b- With gauge conditions K = and D b ^/ ab — 0, as before, we have 

Ah ab = 2 (s ab - p ab ^ cd S cd ] - \labD e h cd D e h cd , (48) 

where £ a t is given by 

Sab ■= -£U^ 4 Aab)+Rat+ 3 Rib~- D aDba-2ip 4 A ac A b c ~8TrS ab .(W) 
a a 

In Ea. H49(l . the first term appears instead of the last three lines in Eq. (|46|) . 

We find that the helically symmetric code does not converge in this method 

when we solve the above set of equations on £ with a boundary that extends several 

wavelengths (or more) beyond the source. We were, however, able to obtain a 

converged solution when exact helical symmetry is imposed only in the near zone, 

within about a wavelength from the source, and the waveless approximation is used 

for larger r, effectively setting boundary conditions at the boundary of the helically 

symmetric inner zone. In a waveless approximation the time derivative of the 

conformal metric, dt^ab is assumed to be zero. As a result the extrinsic curvature is 

associated with the nonrotating shift [3 a , 

K ab = --£„7afc = — £/37a&> (50) 
I la 

instead of the rotating shift ui a as in (|35|l . 

In the next section, we compare the near-zone helically symmetric solution to a 

solution that is everywhere waveless. For the near-zone helically symmetric solution, 

the change from helical symmetry to the waveless formulation at about one wavelength 

from the source implies for K ab the condition 

11 7T 
^-£^7q6 = 7T (£/37ah + fi£tf>7af>) , for r </o> 

la la \l 

K ab = { (51) 

11 7T 

^£/37afc + —labD c (Q,(j) c ), for r > /-, 

where tt/Q, is the approximate wavelength of the I = m = 2 gravitational wave mode. 
The constant /, the coordinate radius of the helically symmetric zone in units of w/Q, 
is restricted to / <, 1 for convergence. 



6.4- Numerical code 

Our numerical code is based on the finite difference code developed in Rcfs. |12l HOj . 
The code extends a KEH iteration scheme to the binary neutron star computation. 
Cartesian components of the field equations are solved numerically on spherical 
coordinate grids, r, 6, and <f>. An equally spaced grid is used from the center of 
orbital motion to 5Rq where there are 16 or 24 grid points per Rq; from 5i?o to the 
outer boundary of computational region a logarithmically spaced grid has 60 or 90, 
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where i?o is the coordinate radius of a compact star along a line passing through the 
center of orbit to the center of a star. Accordingly, for 9 and 4> there are 32 or 48 grid 
points each from to it/ 2 and multipoles are summed up to I — 32 |12| . 

6.5. Numerical solution 

In this section we present results of our code for a binary system modeled by a 
perfect fluid having polytropic equation of state, p = Kp T , with p the baryon density. 
We display results for the choices T = 2, appropriate to neuton star matter; for 
compactness of a star in isolation (M/R)^ — 0.14; and for half the binary separation 
d/Ra = 1.375. Solutions with helical symmetry are not uniquely specified by this 
choice of parameters and equation of state. Because they are stationary solutions 
with equal amounts of ingoing and outgoing waves, they are not asymptotically flat, 
and one must find an appropriate choice of boundary conditions. As discussed in the 
Consortium papers, one seeks conditions that minimize the amplitude of gravitational 
waves. In toy models discussed above, conditions are fixed by the choice of a half- 
advanced plus half-retarded Green function at each iteration. In solving the Einstein 
equation, however, convergence is achieved only in the near zone, and, as we have 
emphasized, we impose boundary conditions by matching to an exterior waveless 
solution outside a coordinate radius /7r/£l. 

For / <, 1, the code converges, yielding a helically symmetric solution in the near 
zone r < fir/fl. As shown in Fig. [7| the solution is nearly identical to the waveless 
solution. The right panel shows a difference larger than 1% only when the metric 
component itself is smaller than 0.03; as a percentage of hij(r = 0), the difference 
is everywhere less than 1%. The threshold of the value of / for convergence is 
0.7 <, f <J 1 depending on the binary separation, compactness and resolution of finite 
differencing. We may expect from the result that, with boundary conditions that 
minimize the amplitude of gravitational waves, the exact helical solution will be close 
to the waveless solution near the source. This is a hopeful outcome: The waveless and 
helically symmetric formalisms are each intended to give a solution whose inaccuracy 
arises from neglecting gravitational waves, and they should give nearly identical results 
in the near zone where the gravitational wave amplitude is small compared to the 
Coulomb fields. 

In modeling binary neutron stars comparable accuracy is likely to result from 
codes that match a helically symmetric solution to a waveless solution, from a purely 
waveless code, and from a helically symmetric code. From a more mathematical 
perspective, however, finding a solution that has exact helical symmetry on the full 
spacetime is an appealing goal. In the neutron star code, we have isolated a set 
of nonlinear terms that appear responsible for divergence outside the near zone. 
Improvement of the convergence of the code is required to extend the matching radius 
beyond a few wavelengths, which may involve a use of metric components in spherical 
instead of Cartesian coordinates. Further investigation of this alternative is a subject 
of our future work. 
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Figure 7. Left panel: Plot of components hij along the x-axis normalized by 
7r/Q, the wavelength of the I = m = 2 mode. The solution to Eq. 1481 with the 
mixed helical and waveless source of Eqs. 1491 and 1511 . and the h yy component 
of the waveless solution are shown. Right panel: The fractional difference of 
these two solutions is plotted for selected components of hij for r < 107r/f2. The 
fractional differences increase for larger r because the values of the components 
hij are themselves small. A compact star extends from x/(tt/Q.) = 0.0125 to 
0.079, the boundary of computational region is set to 10 4 R ~ 332vr/a, and the 
cutoff constant / in Eq. 1511 is set to / = 0.7. 
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